This is the documentation of the \texttt{MODULE Statistics}, a set
of \texttt{FORTRAN 90} routines to perform statistical description
of data. This module make use of the \texttt{MODULE NumTypes},
\texttt{MODULE Constants}, \texttt{MODULE Error} and \texttt{MODULE
  Linear} so please read the documentation of these modules
\emph{before} reading this. 

\section{Function \texttt{Mean(X)}}
\index{Mean@Function \texttt{Mean(X)}}

\subsection{Description}

Compute the mean value of the numbers stored in \texttt{X(:)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose mean we want to compute.
\end{description}

\subsection{Output}

A real double or simple precision (same type as the input). The mean
of the values.

\subsection{Examples}

\begin{lstlisting}[emph=Mean,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Mean of a vector of numbers.,
                   label=mean]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100
  Real (kind=DP) :: X(Nmax)

  CALL Random_Number(X)
  Write(*,'(ES33.25)')Mean(X)

  Stop
End Program Tests
\end{lstlisting}

\section{Function \texttt{Median(X)}}
\index{Median@Function \texttt{Median(X)}}

\subsection{Description}

Compute the median value of the numbers stored in \texttt{X(:)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose median we want to compute.
\end{description}

\subsection{Output}

A real double or simple precision (same type as the input). The median
of the values.

\subsection{Examples}

\begin{lstlisting}[emph=Median,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Median of a vector of numbers.,
                   label=mean]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 5
  Real (kind=SP) :: X(Nmax) = (/1.0, 1.0, 2.0, 4.0, 1.5/)  

  Write(*,'(ES33.25)')Median(X)


  Stop
End Program Tests
\end{lstlisting}

\section{Function \texttt{WMedian(X, w)}}
\index{WMedian@Function \texttt{WMedian(X, w)}}

\subsection{Description}

Compute the weighted median of the numbers stored in \texttt{X(:)}
with weights \texttt{w(:)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose median we want to compute.
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The weights.
\end{description}

\subsection{Output}

A real double or simple precision (same type as the input). The
weighted median
of the values.

\subsection{Examples}

\begin{lstlisting}[emph=WMedian,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Weighted Median of a vector of numbers.,
                   label=wmean]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 5
  Real (kind=SP) :: X(Nmax) = (/1.0, 1.0, 2.0, 4.0, 1.5/)  
  Real (kind=SP) :: w(Nmax) = (/10.0, 2.0, 3.0, 4.0, 1.5/)  

  Write(*,'(ES33.25)')WMedian(X, w)


  Stop
End Program Tests
\end{lstlisting}

\section{Function \texttt{WPercentile(X, w, p)}}
\index{WPercentile@Function \texttt{WPercentile(X, w, p)}}

\subsection{Description}

Compute the weighted percentile \texttt{p} of the numbers stored in
\texttt{X(:)} 
with weights \texttt{w(:)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose median we want to compute.
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The weights.
\item[\texttt{p}:] Double (DP) or simple (SP) precision number. The
  percentile (should be between 0 and 100).
\end{description}

\subsection{Output}

A real double or simple precision (same type as the input). The
weighted percentile \texttt{p}
of the values.

\subsection{Examples}

\begin{lstlisting}[emph=WPercentile,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Weighted Median in two ways.,
                   label=wpercentile]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 5
  Real (kind=SP) :: X(Nmax) = (/1.0, 1.0, 2.0, 4.0, 1.5/)  
  Real (kind=SP) :: w(Nmax) = (/10.0, 2.0, 3.0, 4.0, 1.5/)  

  Write(*,'(ES33.25)')WMedian(X, w)
  Write(*,'(ES33.25)')WPercentile(X, w, 50.0_DP)


  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{WConfInt(X, w, Xmin, Xmax)}}
\index{WConfInt@Subroutine \texttt{WConfInt(X, w, Xmin, Xmax)}}

\subsection{Description}

Compute the weighted 16$^{th}$ and 84$^{th}$ percentiles f the numbers stored in
\texttt{X(:)} 
with weights \texttt{w(:)}, and return the values in \texttt{Xmin} and
\texttt{Xmax}. 

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose median we want to compute.
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The weights.
\item[\texttt{Xmin}:] Double (DP) or simple (SP) precision
  number. Output. The
  16$^{th}$ percentile.
\item[\texttt{Xmax}:] Double (DP) or simple (SP) precision
  number. Output. The
  84$^{th}$ percentile.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=WConfInt,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the 1 sigma confidence interval in two ways., 
                   label=wconfint]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 500
  Real (kind=SP) :: X(Nmax), w(Nmax), X1, X2

  CALL Normal(X)
  CALL Random_Number(w)

  Write(*,'(ES33.25)')WPercentile(X, w, 16.0_DP)
  Write(*,'(ES33.25)')WPercentile(X, w, 84.0_DP)
  CALL WConfInt(X, W, X1, X2)
  Write(*,'(ES33.25)')X1, X2
  
  Stop
End Program Tests
\end{lstlisting}

\section{Function \texttt{Var(X)}}
\index{Var@Function \texttt{Var(X)}}

\subsection{Description}

Compute the variance of a vector of numbers \texttt{X(:)}

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose variance we want to compute.
\end{description}

\subsection{Output}

A real double or simple precision (same type as the input). The
variance of the values.

\subsection{Examples}

\begin{lstlisting}[emph=Var,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Variance of a set of numbers.,
                   label=var]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100, Npinta = 100, Npar = 4
  Real (kind=DP) :: X(Nmax), Y(Nmax), Yer(Nmax), &
       & Coef(Npar), Cerr(Npar), Corr, Xd(Nmax,2)


  CALL Random_Number(X)
  Write(*,'(ES33.25)')Var(X)


  Stop
End Program Tests
\end{lstlisting}


\section{Function \texttt{Stddev(X)}}
\index{Stddev@Function \texttt{Stddev(X)}}

\subsection{Description}

Computes the standard deviation of the numbers stored in the vector
\texttt{X(:)}. 

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Double (DP) or simple (SP) precision one
  dimensional array. The values  whose standard deviation we want to
  compute. 
\end{description}

\subsection{Output}

Real Single or Double precision, the same as the input values. The
standard deviation of the values.

\subsection{Examples}

\begin{lstlisting}[emph=Stddev,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Compputing the standard deviation.,
                   label=stddev]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100, Npinta = 100, Npar = 4
  Real (kind=DP) :: X(Nmax), Y(Nmax), Yer(Nmax), &
       & Coef(Npar), Cerr(Npar), Corr, Xd(Nmax,2)


  CALL Random_Number(X)
  Write(*,'(ES33.25)')Stddev(X)


  Stop
End Program Tests
\end{lstlisting}

\section{Function \texttt{Moment(X, k)}}
\index{Moment@Function \texttt{Moment(X, k)}}

\subsection{Description}

Returns the $k^{\underline{th}}$ moment of the values stored in the
vector \texttt{X(:)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision). The numbers
  whose $k^{\underline{th}}$ moment we want to compute.
\item[\texttt{k}:] Integer. Which moment we want to compute.
\end{description}

\subsection{Output}

Real single or double precision. The $k^{\underline{th}}$ moment of
the numbers.

\subsection{Examples}

\begin{lstlisting}[emph=Moment,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the k$^{\text{\underline{th}}}$
                   moment of a data set.,
                   label=moment]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100, Npinta = 100, Npar = 4
  Real (kind=DP) :: X(Nmax), Y(Nmax), Yer(Nmax), &
       & Coef(Npar), Cerr(Npar), Corr, Xd(Nmax,2)


  CALL Random_Number(X)
  Write(*,*)'We should obtain the same numbers twice: '
  Write(*,'(ES33.25)')Moment(X,2), Var(X)

  Stop
End Program Test
\end{lstlisting}

\section{Subroutine \texttt{Histogram(Val, Ndiv, Ntics, Vmin, Vmax,
    h)}} 
\index{Histogram@Subroutine \texttt{Histogram(Val, Ndiv, Ntics, Vmin, Vmax, h)}} 

\subsection{Description}

Given a set of points \texttt{Val(:)}, this routine makes
\texttt{Ndiv} divisions between the minimum and the greatest value of
\texttt{Val} (respectively returned in \texttt{Vmin} and
\texttt{Vmax}), each of size \texttt{h} (also returned), and returns
in the integer vector \texttt{Nticks(:)} the number of points that are
in each interval. 

\subsection{Arguments}

\begin{description}
\item[\texttt{Val(:)}:] Real (Single or Double precision) one
  dimensional array. The original values.
\item[\texttt{Ndiv}: ] Integer. The number of divisions.
\item[\texttt{Nticks}:] Integer one dimensional array. \texttt{Ndiv(I)}
  Tells how many points of \texttt{Val(:)} are between
  $\mathtt{Vmin+(I-1)h}$ and $\mathtt{Vmin+Ih}$.
\item[\texttt{Vmin, Vmax}:] Real (Single or Double precision). The
  minimum and maximum values of \texttt{Val}.
\item[\texttt{h}:] Real (Single or Double precision). After calling
  the routine has the step of the division.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Histogram,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Making Histograms.,
                   label=histogram]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 500000, Npinta = 100, Npar = 4, Ndiv = 100
  Real (kind=DP) :: X(Nmax), Y(Nmax), Yer(Nmax), &
       & Coef(Npar), Cerr(Npar), Corr, Xd(Nmax,2), &
       & Xmin, Xmax, h, Xac
  Integer :: Ntics(Ndiv)

  CALL Normal(X, 1.23_DP, 0.345_DP)
  CALL Histogram(X, Ndiv, Ntics, Xmin, Xmax, h)
  
  Do I = 1, Ndiv
     Xac = Xmin + (I-1)*h
     Write(*,'(1ES33.25,1I)')Xac, Ntics(I)
  End Do

  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{LinearReg(X, Y, Yerr, [Func], Coef, Cerr, ChisqrV)}} 
\index{LinearReg@Subroutine \texttt{LinearReg(X, Y, Yerr, [Func], Coef, Cerr, ChisqrV)}} 

\subsection{Description}

Given a set of points \texttt{X(:)} and \texttt{Y(:)}, this routine
performs a linear fit to a set of functions defined by
\texttt{Func}. 
\begin{displaymath}
  Y = \sum_i a_i f_i(X)
\end{displaymath}
This routine also performs multi-dimensional fitting, in which case
the points are specified as \texttt{X(:,:)}, where the first argument
tells which point, and the second which variable.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:[,:])}:] Real single or double precision one
  dimensional array (for a one dimensional fit) or two dimensional
  array (for a multidimensional fit). The
  independent variables. For a multidimensional fit, the first argument
  tells which point, and the second which variable. So the size of the
  array should be \texttt{X(Npoints,Ndim)}.
\item[\texttt{Y(:)}: ] Real single or double precision one dimensional
  array. The dependent
  variable.
\item[\texttt{Yerr(:[,:])}:] Real single or double precision one or
  two
  dimensional array. If a one dimensional array is inserted, they are
  the errors 
  of the points (if you don't have them, you should put all of them to
  some non-zero value). If a two dimensional array is given, it is
  treated as the correlation matrix.
\item[\texttt{Func}:] Optional. This routine define the functions to
  fit. An interface like this should be provided
\begin{verbatim}
Interface
   Function Func(Xx, i)
         
     USE NumTypes

     Real (kind=SP), Intent (in) :: Xx
     Integer, Intent (in) :: i
     Real (kind=SP) :: Func

   End Function Func
End Interface
\end{verbatim}
if you want to perform a one dimensional fitting, and like this
\begin{verbatim}
Interface
   Function Func(Xx, i)
         
     USE NumTypes

     Real (kind=SP), Intent (in) :: Xx(:)
     Integer, Intent (in) :: i
     Real (kind=SP) :: Func

   End Function Func
End Interface
\end{verbatim}
if it is a multidimensional fitting. Since you are making a fitting
to a function of the type
\begin{displaymath}
  Y = \sum_i a_i f_i(X)
\end{displaymath}
the values $f_i(X)$ are given by this function as \texttt{Func(X,
  I)}. If the functions are not specified (i.e. you don't put this
argument), a fit to a polynomial is made (this only work for
one-dimensional fittings).
\item[\texttt{Coef(:)}: ] Real single or double precision one
  dimensional array. The
  parameters that you want to determine.
\item[\texttt{Cerr(:)}:] Real single or double precision one
  dimensional array. The errors
  in the parameters.
\item[\texttt{ChiSqr}: ] Real single or double precision. The $\chi^2$
  per degree of freedom of the fit.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=LinearReg,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Doing linear regressions.,
                   label=linearreg]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 200, Npinta = 100, Npar = 4, Ndiv = 100
  Real (kind=DP) :: X(Nmax), Y(Nmax), Yer(Nmax), &
       & Coef(Npar), Cerr(Npar), Corr, Xd(Nmax,2), &
       & Xmin, Xmax, h, Xac
  Integer :: Ntics(Ndiv)

  Interface
     Function Fd(Xx, i)
       
       USE NumTypes
       
       Real (kind=DP), Intent (in) :: Xx(:)
       Integer, Intent (in) :: i
       Real (kind=DP) :: Fd
       
     End Function FD
  End Interface


  CALL Random_Number(Xd)
  Xd(:,:) = 10.0_DP*(Xd(:,:) - 0.8_DP)

  CALL Normal(Yer, 0.0_DP, 1.0E-3_DP)
  Y(:) = 12.34_DP*Xd(:,1)*sin(Xd(:,2)) - 2.23_DP + &
       & 0.67_DP*Xd(:,1)**2*Xd(:,2) +  0.23_DP*Xd(:,1) + Yer(:) 


  CALL LinearReg(Xd, Y, Yer, Fd, Coef, Cerr, Corr)
  
  ! This should print the adjusted parameters, 
  ! that have values: 12.34, -2.23, 0.67, 0.23
  Do I = 1, Npar
     Write(*,'(2ES33.25)')Coef(I), Cerr(I)
  End Do

  ! This prints the ChiSqr, that should be very 
  ! close to 1.
  Write(*,'(1A,1ES33.25)')'ChiSqr of the Fit: ', Corr


  Stop
End Program Tests

! ************************************
! *
Function Fd(X, i)
! *
! ************************************

  USE NumTypes

  Real (kind=DP), Intent (in) :: X(:)
  Integer, Intent (in) :: i
  Real (kind=DP) :: Fd

  If (I==1) Then
     Fd = 1.0_DP
  Else If (I==2) Then
     Fd = X(1)*sin(X(2))
  Else If (I==3) Then
     Fd = X(1)**2*X(2)
  Else If (I==4) Then
     Fd = X(1)
  End If

  Return
End Function FD
\end{lstlisting}


\section{Subroutine \texttt{NonLinearReg(X, Y, Yerr, Func, Coef, Cerr, ChisqrV)}} 
\index{NonLinearReg@Subroutine \texttt{NonLinearReg(X, Y, Yerr, Func, Coef, Cerr, ChisqrV)}} 

\subsection{Description}

Given a set of points \texttt{X(:)} and \texttt{Y(:)}, this routine
performs a non-linear fit to a set of functions defined by
\texttt{Func}.  

This routine also performs multi-dimensional fitting, in which case
the points are specified as \texttt{X(:,:)}, where the first argument
tells which point, and the second which variable.

This routine uses the Levenberg-Marquardt algorithm to perform the
optimisation\footnote{\href{http://en.wikipedia.org/wiki/Levenberg-Marquardt_algorithm}{\texttt{http://en.wikipedia.org/wiki/Levenberg-Marquardt\_algorithm}}}. 

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:[,:])}:] Real single or double precision one
  dimensional array (for a one dimensional fit) or two dimensional
  array (for a multidimensional fit). The
  independent variables. For a multidimensional fit, the first argument
  tells which point, and the second which variable. So the size of the
  array should be \texttt{X(Npoints,Ndim)}.
\item[\texttt{Y(:)}: ] Real single or double precision one dimensional
  array. The dependent
  variable.
\item[\texttt{Yerr(:[,:])}:] Real single or double precision one or
  two
  dimensional array. If a one dimensional array is inserted, they are
  the errors 
  of the points (if you don't have them, you should put all of them to
  some non-zero value). If a two dimensional array is given, it is
  treated as the correlation matrix.
\item[\texttt{Func}:] This routine define the functions to
  fit. An interface like this should be provided
\begin{verbatim}
Interface
   Subroutine Func(X, Cf, Valf, ValD)
         
     USE NumTypes

     Real (kind=SP), Intent (in) :: X, Cf(:)
     Real (kind=SP), Intent (out) :: Valf, ValD(Size(Cf))
         
   End Subroutine Func
End Interface
\end{verbatim}
if you want to perform a one dimensional fitting, and like this
\begin{verbatim}
Interface
   Subroutine Func(X, Cf, Valf, ValD)
        
     USE NumTypes

     Real (kind=SP), Intent (in) :: X(:), Cf(:)
     Real (kind=SP), Intent (out) :: Valf, ValD(Size(Cf))
         
   End Subroutine Func
End Interface
\end{verbatim}
if it is a multidimensional fitting. 

This routine returns the values of the function at $X$ for some values
of the parameters given in $\mathtt{Cf}$ in the variable
$\mathtt{Valf}$, and a vector with the derivatives (respect with the
parameters) in $\mathtt{ValD(:)}$. 
\item[\texttt{Coef(:)}: ] Real single or double precision one
  dimensional array. Output. The
  parameters that you want to determine.
\item[\texttt{Cerr(:)}:] Real single or double precision one
  dimensional array. Outpue. The errors
  in the parameters.
\item[\texttt{ChiSqr}: ] Real single or double precision. The $\chi^2$
  per degree of freedom of the fit.
\end{description}

\subsection{Examples}

In this example we will fit some generated data, with a Normal noise
to the function
\begin{displaymath}
  f(x_1,x_2;a,b) = \sin(ax_1) + bx_1x_2
\end{displaymath}
We will generate the data with the values $a=2.0$ and $b=0.2$, so our
fitting routine \emph{should} return this values within errors.

The derivatives of the function (repect the parameters $a$ and $b$),
as well as the value of function are given by the user routine
\texttt{Func}. The derivatives are:
\begin{eqnarray*}
  \frac{\partial f(x_1,x_2;a,b)}{\partial a} &=& x_1\cos{ax_1} \\
  \frac{\partial f(x_1,x_2;a,b)}{\partial b} &=& x_1x_2 \\
\end{eqnarray*}


\begin{lstlisting}[emph=NonLinearReg,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Doing non linear regressions.,
                   label=linearreg]
Program NLFit

  USE NumTypes
  USE Statistics

  Interface
     Subroutine Func(X, Cf, Valf, ValD)
       
       USE NumTypes
       
       Real (kind=DP), Intent (in) :: X(:), Cf(:)
       Real (kind=DP), Intent (out) :: Valf, ValD(Size(Cf))
       
     End Subroutine Func
  End Interface

  Integer, Parameter :: Np = 20, Ndim = 2
  Real (kind=DP) :: X(Np, Ndim), Y(Np), Ye(Np), Co(2), Vd(2), Ce(2), Ch

  ! First Fill the data
  CALL Random_Number(X)
  X = 2.0_DP*(X - 0.5_DP)
  Co(1) = 2.0_DP
  Co(2) = 0.2_DP
  CALL Normal(Ye, 0.0_DP, 0.5_DP)
  Do I = 1, Np
     CALL Func(X(I,:), Co, Y(I), Vd)
     Y(I) = Y(I) + Ye(I)
  End Do
  

  ! Now Perform the non linear fit
  Co = 1.0_DP
  CALL NonLinearReg(X, Y, Abs(Ye), Func, Co, Ce, Ch)
  Do I = 1, Npar
    Write(*,'(1A,100ES33.25)')'Parameter and error: ', Co(I), Ce(I)
  End Do
  Write(*,'(1A,100ES33.25)')&
  & 'Chi Square per degree of freedom of the Fit: ', Ch


  Stop
End Program NLFit

Subroutine Func(X, Cf, Valf, ValD)
  
  USE NumTypes
  
  Real (kind=DP), Intent (in) :: X(:), Cf(:)
  Real (kind=DP), Intent (out) :: Valf, ValD(Size(Cf))

  Valf = Sin(Cf(1)*X(1)) + Cf(2)*X(1)*X(2)
  ValD(1) = X(1)*Cos(Cf(1)*X(1))
  ValD(2) = X(1)*X(2)

  
End Subroutine Func
\end{lstlisting}

\section{Subroutine \texttt{SetLuxLevel(Ilevel)}}
\index{SetLuxLevel@Subroutine \texttt{SetLuxLevel(Ilevel)}}

\subsection{Description}

This function changes the behaviour of the ``intrinsic''
\texttt{Random\_Number} subroutine. To call the intrinsic
\texttt{FORTRAN 90} routine, you should set the Luxury level to
zero. Other values makes the \texttt{Random\_Number} subroutine use the
Marsaglia and Zaman algorithm modified by
M. L\"uscher~\cite{Luscher:1993dy}. This part of the code has some
inspiration in the  \texttt{FORTRAN 90} implementation of Allan
Miller. \textbf{Note:} This routine affects the way all the other
routines that use pseudo random number generators works. That is to
say, all the routines of this module that follow this one.

\subsection{Arguments}

\begin{description}
\item[\texttt{Ilevel}:] Integer. Set the way the random number
  generator works. We have the following options:
  \begin{itemize}
  \item \texttt{Ilevel} $=0$. The default intrinsic \texttt{FROTRAN
      90} routine is used. This is usually considered a \emph{bad},
    non portable pseudo-random number generator, that you do not want
    to use when correlations are important.
  \item \texttt{Ilevel} $=1$. This is the original Marsaglia and
    Zaman algorithm. Probably better than the intrinsic procedure, but
    still with large correlations.
  \item \texttt{Ilevel} $=2$. The M. L\"uscher modified version of the
    algorithm with $p=48$. Still fail many tests, but a very good
    pseudo-number generator for many things.
  \item \texttt{Ilevel} $=3$. The M. L\"uscher modified version of the
    algorithm with $p=97$. Theoretically still defective (fail the
    serial correlation test).
  \item \texttt{Ilevel} $=4$. The M. L\"uscher modified version of the
    algorithm with $p=218$. Any theoretically possible correlation has
    a very small probability of being observed. A value very similar
    to this has been
    used for large scale high precision lattice gauge theory
    computations. This corresponds with the L\"uscher recommended level
    0. 
  \item \texttt{Ilevel} $=5$. The M. L\"uscher modified version of the
    algorithm with $p=404$. Good for
    anything you need the pseudo-random numbers (but cryptography, of
    course). It seems completely pointless to use values of $p$ higher
    than this.  This corresponds with the L\"uscher recommended level
    1. This is the \textbf{default value}.
  \item \texttt{Ilevel} $=6$. The M. L\"uscher modified version of the
    algorithm with $p=794$. Good for
    anything you need the pseudo-random numbers (but cryptography, of
    course). Probably waste of computer resources.  This corresponds
    with the L\"uscher recommended level 2. 
  \item Other values of \texttt{Ilevel} $(>24)$. Sets the value of $p$
    in the M. L\"uscher modified version of the algorithm. Don't use
    it, unless you know what you are doing. 
  \end{itemize}
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=SetLuxLevel,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Setting the Luxury level of the pseudo random number generator., 
                   label=setluxlevel]
Program Tests

  USE NumTypes 
  USE Statistics

  Integer, Parameter :: NN = 70
  Real (kind=SP) :: rr(NN)

  CALL SetLuxLevel(5)
  CALL Random_Number(rr)
  Write(*,*)rr

  
  Stop
End Program Tests
\end{lstlisting}


\section{Subroutine \texttt{PutLuxSeed([ISeed(25)])}}
\index{PutLuxSeed@Subroutine \texttt{PutLuxSeed([ISeed])}}

\subsection{Description}

Restarts the position in the Luxury pseudo-random number generator. 

\subsection{Arguments}

\begin{description}
\item[\texttt{ISeed(25)}:] Integer one dimensional array of 25
  elements. Optional. If present, restarts the generator from an output
  of \texttt{GetLuxSeed}. If not present, initialises the random
  number generator.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=PutLuxSeed,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using a previously saved point in the generating process.,
                   label=putluxseed]
Program Tests

  USE NumTypes 
  USE Statistics

  Integer, Parameter :: NN = 70
  Real (kind=SP) :: rr(NN)
  Integer :: Sd(25)


  Open(Unit=69, File="seed.dat")
  Read(69,*)Sd
  CALL PutLuxSeed(Sd)
  CALL Random_Number(rr)
  Write(*,*)rr

  
  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{GetLuxSeed(ISeed(25))}}
\index{GetLuxSeed@Subroutine \texttt{GetLuxSeed([ISeed])}}

\subsection{Description}

Saves the position in the Luxury pseudo-random number generator. 

\subsection{Arguments}

\begin{description}
\item[\texttt{ISeed(25)}:] Integer one dimensional array of 25
  elements. Saves the position in the random number generator.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=GetLuxSeed,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Saving a point in the generating process.,
                   label=getluxseed]
Program Tests

  USE NumTypes 
  USE Statistics

  Integer, Parameter :: NN = 70
  Real (kind=SP) :: rr(NN)
  Integer :: Sd(25)


  CALL Random_Number(rr)
  CALL GetLuxSeed(Sd)
  Open(Unit=69, File="seed.dat")
  Write(69,*)Sd

  
  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{WriteLuxSeed(fname,[comment[(:)]])}}
\index{WriteLuxSeed@Subroutine \texttt{WriteLuxSeed(fname,[comment[(:)]])}}

\subsection{Description}

Writes the position of the LUX random number generator in file
\texttt{fname}. 

\subsection{Arguments}

\begin{description}
\item[\texttt{fname}:] Character (len=*). The file name.
\item[\texttt{Comment[(:)]}:] Character (len=*) or array of Character
  (len=*). Comments to add to the header of the file. 
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=WriteLuxSeed,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using a previously saved point in the generating process.,
                   label=writeluxseed]

Program Seed

  USE NumTypes
  USE Statistics

  Real (kind=DP), Allocatable :: X(:)
  Character (len=50) :: CC(5)

  CC(1) = 'Follow Input'
  CC(2) = 'a=1'
  CC(3) = 'b=2'
  CC(3) = 'c=23'
  CC(3) = 'd=276278'

  CALL PutLuxSeed()
  Allocate(X(1000))
  CALL Normal(X)
  CALL WriteLuxSeed('o.seed', CC)

  CALL Normal(X)
  Write(*,*)X(990:1000)

  Write(*,*)'Again: '
  CALL ReadLuxSeed('o.seed')
  CALL Normal(X)
  Write(*,*)X(990:1000)

  Stop
End Program Seed
\end{lstlisting}


\section{Subroutine \texttt{ReadLuxSeed(fname)}}
\index{ReadLuxSeed@Subroutine \texttt{ReadLuxSeed(fname)}}

\subsection{Description}

Reads the position of the LUX random number generator from file
\texttt{fname}.

\subsection{Arguments}

\begin{description}
\item[\texttt{fname}:] Character (len=*). The file name.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=ReadLuxSeed,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using a previously saved point in the generating process.,
                   label=readluxseed]

Program Seed

  USE NumTypes
  USE Statistics

  Real (kind=DP), Allocatable :: X(:)
  Character (len=50) :: CC(5)

  CC(1) = 'Follow Input'
  CC(2) = 'a=1'
  CC(3) = 'b=2'
  CC(3) = 'c=23'
  CC(3) = 'd=276278'

  CALL PutLuxSeed()
  Allocate(X(1000))
  CALL Normal(X)
  CALL WriteLuxSeed('o.seed', CC)

  CALL Normal(X)
  Write(*,*)X(990:1000)

  Write(*,*)'Again: '
  CALL ReadLuxSeed('o.seed')
  CALL Normal(X)
  Write(*,*)X(990:1000)

  Stop
End Program Seed
\end{lstlisting}

\section{Function \texttt{LUXRandomize([N])}}
\index{LUXRandomize@Function \texttt{LUXRandomize(N)}}

\subsection{Description}

Randomize the LUX random number generator based on actual time or a
reduced seed. 

\subsection{Arguments}

\begin{description}
\item[\texttt{N}:] Integer, Optional. If present a reduced seed. If
  not present it is generated based on hour, minutes, seconds and
  miliseconds of the calling time.
\end{description}

\subsection{Output}

Integer. The ``reduced seed'' used. 

\subsection{Examples}

\begin{lstlisting}[emph=LUXRandomize,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Randomizing the LUX generator., 
                   label=luxrandomize]

Program Seed

  USE NumTypes
  USE Statistics

  Real (kind=DP), Allocatable :: X(:)

  CALL PutLuxSeed()
  Iseed = LUXRandomize()
  Allocate(X(1000))
  CALL Normal(X)

  Stop
End Program Seed
\end{lstlisting}


\section{Subroutine \texttt{Normal(X, [Rm], [Rsig])}}
\index{Normal@Subroutine \texttt{Normal(X, [Rm], [Rsig])}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from a normal distribution with mean
\texttt{Rm}, and standard deviation \texttt{Rsig}. The parameters
\texttt{Rm} and \texttt{Rsig} are optional. If they are not given the
mean will be 0, and the standard deviation 1.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the normal distribution.
\item[\texttt{Rm}:] Real (Single or Double precision), Optional. The
  mean of the normal distribution. If not present the default value
  is 0.
\item[\texttt{Rsig}:] Real (Single or Double precision), Optional. The
  standard deviation of the normal distribution.  If not present the
  default value is 1. 
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Normal,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining numbers with a normal distribution.,
                   label=normal]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100
  Real (kind=DP) :: X(Nmax)


  CALL Normal(X, 1.23_DP, 0.345_DP)
  ! Now compute the mean and standard deviation of the data
  Write(*,*)'We should obtain 1.23 and 0.345: '
  Write(*,'(ES33.25)')Mean(X), Stddev(X)


  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{MultiNormal(X, [Rm], Sig)}}
\index{MultiNormal@Subroutine \texttt{MultiNormal(X, [Rm], Sig)}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from a multivariate normal
distribution with mean (optional) \texttt{Rm(:)}, and covariance matrix
\texttt{Sig(:)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:[,:])}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the multivariate normal distribution with covariance
  matrix \texttt{Sig}. If a two dimensional array is introduced of
  dimensions $N\times N_v$ the array will be filled with $N$ samples
  of $N_v$ multivariate normally distributed numbers with covariance
  given by the $N_v\times N_v$ matrix \texttt{Sig(:,:)}.
\item[\texttt{Rm(:)}:] Real (Single or Double precision), Optional. The
  means of the multivariate vector. If not present the default value
  is 0.
\item[\texttt{Rsig(:,:)}:] Real (Single or Double precision). The
  Covariance matrix. 
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=MultiNormal,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Sampling a multivariate normal distribution.,
                   label=multinormal]
Program MultiN

  USE NumTypes
  USE Statistics

  Integer, Parameter :: N = 2, Ns = 20000000
  Real (kind=DP) :: X(N), Z(Ns, N), S(N,N), r, A(N,N)

  S(1,1) = 4.0_DP
  S(1,2) = 0.2_DP
  S(2,1) = 0.2_DP
  S(2,2) = 1.0_DP


  ! Z(:,1) and Z(:,2) are normally sampled with covariance S
  CALL MultiNormal(Z, (/1.0_DP, 2.0_DP/), S)

  r = 0.0_DP
  Do I = 1, Ns
     r = r + (Z(I,1)-1.0_DP)/2.0_DP * (Z(I,2)-2.0_DP)
  End Do
  r = r / Real(Ns-1,kind=DP)
  
  Write(*,*)"Pearson's correlation coefficient: ", r

  Stop
End Program MultiN
\end{lstlisting}



\section{Subroutine \texttt{Cauchy(X, [c], [$\beta$], [$\mu$])}}
\index{Cauchy@Subroutine \texttt{Cauchy(X, [c], [$\beta$], [$\mu$])}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from a stable distribution with
mean $\mu$, scale c and exponent $\alpha=1$. The probability
distribution is defined via
\begin{equation}
  \rho(x) = \frac{1}{\pi} \Re\int_{-\infty}^{+\infty} dt 
  \exp\left\{ -\imath t (x-\mu) - c|t|^\alpha\left[1+\imath\frac{t\beta}{|t|}\tan\left(\frac{\pi\alpha}{2}\right)\right] \right\}
\end{equation}

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the normal distribution.
\item[$\alpha$:] Real (Single or Double precision). The
  characteristic exponent (or Levy index).
\item[\texttt{c}:] Real (Single or Double precision), Optional. The
  scale factor.
\item[$\beta$:] Real (Single or Double precision), Optional. The 
  skewness parameter.
\item[$\mu$:] Real (Single or Double precision), Optional. The 
  translation factor.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Cauchy,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining numbers with a stable distribution.,
                   label=cauchy]
Program CacuchyTest

  USE NumTypes
  USE Statistics
  USE NonNumeric

  Integer, Parameter :: N = 5000000, Nd = 200
  Real (kind=DP) :: X(N), Y(N), V1, V2, h, Div(Nd), alph, b, s, m
  Integer :: Nt1(Nd), Nt2(Nd), Id
  

  Read(*,*)b, s, m
 
  CALL Cauchy(X, b, s, m)
  CALL Normal(Y, 0.0_DP, SR2_DP)

  h = 20.0_DP/Real(Nd,kind=DP)
  Do I = 1, Nd
     Div(I) = -10.0_DP + (I-1)*h
  End Do

  Nt1 = 0
  Nt2 = 0
  Do I = 1, N
     Id = Locate(Div, X(I))
     Nt1(Id) = Nt1(Id) +1
     Id = Locate(Div, Y(I))
     Nt2(Id) = Nt2(Id) +1
  End Do

  Do I = 2, Nd-1
     Write(*,*)-10.0_DP + (I-1)*h, Real(Nt1(I))/Real(N), Real(Nt2(I))/Real(N)
  End Do

  Write(0,*)Mean(X), Stddev(X)
  Write(0,*)Mean(Y), Stddev(Y)


  Stop
End Program CacuchyTest
\end{lstlisting}


\section{Subroutine \texttt{Lorentz(X, [$\mu$], [$\gamma$])}}
\index{Lorentz@Subroutine \texttt{Lorentz(X, [$\mu$], [$\gamma$])}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from Lorentz distribution with
location $\mu$ and $\gamma=1$. The probability 
distribution is defined via
\begin{equation}
  \rho(x) = \frac{1}{\pi} \frac{\gamma}{(x-\mu)^2 + \gamma^2}
\end{equation}

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the normal distribution.
\item[$\mu$:] Real (Single or Double precision), Optional. The 
  translation factor.
\item[$\gamma$:] Real (Single or Double precision), Optional. The 
  scale factor.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Lorentz,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Comparing Lorentz and Cauchy routines.,
                   label=lorentz]
Program LorentzTest

  USE NumTypes
  USE Statistics
  USE NonNumeric

  Integer, Parameter :: N = 50000, Nd = 200
  Real (kind=DP) :: X(N), Y(N), Z(N), V1, V2, h, Div(Nd), alph, b, s, m
  Integer :: Nt1(Nd), Nt2(Nd), Nt3(Nd), Id
  

  Read(*,*)m, s
 
  CALL Lorentz(X, m, s)
  CALL Cauchy(Y, 0.0_DP, s, m)

  Z = 2*X/(s**2+X**2)

  h = 20.0_DP/Real(Nd,kind=DP)
  Do I = 1, Nd
     Div(I) = -10.0_DP + (I-1)*h
  End Do

  Nt1 = 0
  Nt2 = 0
  Nt3 = 0
  Do I = 1, N
     Id = Locate(Div, X(I))
     Nt1(Id) = Nt1(Id) +1
     Id = Locate(Div, Y(I))
     Nt2(Id) = Nt2(Id) +1
     Id = Locate(Div, Z(I))
     Nt3(Id) = Nt3(Id) +1
  End Do

  Do I = 2, Nd-1
     Write(*,'(100ES23.15)')-10.0_DP + I*h, &
          & Real(Nt1(I))/Real(N)/h, &
          & Real(Nt2(I))/Real(N)/h, &
          & Real(Nt3(I))/Real(N)/h
  End Do



  Write(0,*)Mean(X), Stddev(X)
  Write(0,*)Mean(Y), Stddev(Y)
  Write(0,*)Sum(Nt3(:)), N

  Stop
End Program LorentzTest
\end{lstlisting}


\section{Subroutine \texttt{Levy(X, $\alpha$, [c], [$\beta$], [$\mu$])}}
\index{Levy@Subroutine \texttt{Levy(X, $\alpha$, [c], [$\beta$], [$\mu$])}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from a Levy stable distribution with
mean $\mu$, scale c and exponent $\alpha$. The probability
distribution is defined via
\begin{equation}
  \rho(x) = \frac{1}{\pi} \Re\int_{-\infty}^{+\infty} dt 
  \exp\left\{ -\imath t (x-\mu) - c|t|^\alpha\left[1+\imath\frac{t\beta}{|t|}\tan\left(\frac{\pi\alpha}{2}\right)\right] \right\}
\end{equation}
For $\alpha = 1$ the
distribution reduces to the Cauchy distribution, but this function is
unstable. The cauchy routine should be used.

The algorithm only works for $0 < \alpha \le 2$ and $\alpha \neq 1$. 

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the normal distribution.
\item[$\alpha$:] Real (Single or Double precision). The
  characteristic exponent (or Levy index).
\item[\texttt{c}:] Real (Single or Double precision), Optional. The
  scale factor.
\item[$\beta$:] Real (Single or Double precision), Optional. The 
  skewness parameter.
\item[$\mu$:] Real (Single or Double precision), Optional. The 
  translation factor.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Levy,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining numbers with a Levy skew stable distribution.,
                   label=normal]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100
  Real (kind=DP) :: X(Nmax)


  CALL Levy(X, 1.23_DP)
  ! Now compute the mean and standard deviation of the data
  Write(*,*)'We should obtain 0.0 and something without sense!: '
  Write(*,'(ES33.25)')Mean(X), Stddev(X)


  Stop
End Program Tests
\end{lstlisting}


\section{Subroutine \texttt{FishTipp(X, Rm, Rb)}}
\index{FishTipp@Subroutine \texttt{FishTipp(X, Rm, Rb)}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from a Fisher-Tippet distribution with
parameters\footnote{More info about this distribution in the
  Wikipedia: \href{http://en.wikipedia.org/wiki/Fisher-Tippett_distribution}{\texttt{http://en.wikipedia.org/wiki/Fisher-Tippett\_distribution}}} \texttt{Rm}, and $2\mathtt{Rb}^2$.

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the normal distribution.
\item[\texttt{Rm}:] Real (Single or Double precision). 
\item[\texttt{Rb}:] Real (Single or Double precision). 
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=FishTipp,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining numbers with a Fisher-Tippet distribution.,
                   label=fishtipp]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100
  Real (kind=DP) :: X(Nmax)


  CALL FishTipp(X, 2.00_DP, 1.00_DP)
  ! Now compute the mean and standard deviation of the data
  Write(*,*)'We should obtain 2.57721... and 1.2782...: '
  Write(*,'(ES33.25)')Mean(X), Stddev(X)


  Stop
End Program Tests
\end{lstlisting}


\section{Subroutine \texttt{Laplace(X, Rm, Rb)}}
\index{Laplace@Subroutine \texttt{Laplace(X, Rm, Rb)}}

\subsection{Description}

Fills \texttt{X(:)} with numbers from a Laplace distribution with mean
\texttt{Rm}, and variance $2\mathtt{Rb}^2$. 

\subsection{Arguments}

\begin{description}
\item[\texttt{X(:)}:] Real (Single or Double precision) one
  dimensional array. A vector that will be filled with numbers
  according to the Laplace distribution.
\item[\texttt{Rm}:] Real (Single or Double precision). The
  mean of the Laplace distribution.
\item[\texttt{Rb}:] Real (Single or Double precision). The
  width of the Laplace distribution (i.e. The variance is
  $2\mathtt{Rb}^2$).
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Laplace,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining numbers with a Laplace distribution.,
                   label=laplace]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100
  Real (kind=DP) :: X(Nmax)


  CALL Laplace(X, 1.23_DP, 1.0_DP)
  ! Now compute the mean and standard deviation of the data
  Write(*,*)'We should obtain 1.23 and sqrt(2): '
  Write(*,'(ES33.25)')Mean(X), Stddev(X)


  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine/Function \texttt{Irand([Irnd], N, M)}}
\index{Irand@Subroutine/Function \texttt{Irand([Irnd], N, M)}}

\subsection{Description}

If present, fills \texttt{Irnd(:)} with random integer numbers between
\texttt{N} and \texttt{M} with an uniform distribution. If
\texttt{Irnd(:)} is not present returns a integer random number
between \texttt{N} and \texttt{M}.

\subsection{Arguments}

\begin{description}
\item[\texttt{Irnd(:)}:] Integer, Optional. A vector that will be
  filled with integer numbers according to a uniform distribution.
\item[\texttt{N}:] Integer. The minimum number that we can obtain.
\item[\texttt{M}:] Integer. The maximum number that we can obtain.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Irand,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining integer random numbers.,
                   label=irand]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100
  Integer :: Irnd(Nmax)


  CALL Irand(Irnd, 0, 1)
  ! Now compute the mean 
  Write(*,*)'We should obtain 0.5: '
  Write(*,'(ES33.25)')Mean(Real(Irnd(:),kind=DP))


  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{Permutation(Idx)}}
\index{Permutation@Subroutine \texttt{Permutation(Idx)}}

\subsection{Description}

Returns a random permutation of $\mathtt{N}$ elements. It uses the
Knuth shuffle algorithm\footnote{\href{http://en.wikipedia.org/wiki/Knuth_shuffle}{\texttt{http://en.wikipedia.org/wiki/Knuth\_shuffle}}}.

\subsection{Arguments}

\begin{description}
\item[\texttt{Idx(:): }] Integer  one dimensional array. Output. The
  random permutation.
\end{description}

\begin{lstlisting}[emph=Permutation,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Obtaining a permutation.,
                   label=Permutation]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 10
  Integer :: Id(Nmax)


  CALL Permutation(Id)
  Write(*,'(100I3)')(Id(I), I = 1, Nmax)

  Stop

End Program Tests
\end{lstlisting}


\section{Subroutine \texttt{BootStrap(Ibt)}}
\index{Bootstrap@Subroutine \texttt{BootStrap(Ibt)}}

\subsection{Description}

Generates $\mathtt{N_b}$ bootstrap sequence of $\mathtt{N}$ numbers
each. These bootstraps are returned in the two dimensional integer
array $\mathtt{Ibt(:,:)}$ of size $\mathtt{N\times N_b}$.

\subsection{Arguments}

\begin{description}
\item[\texttt{Ibt(:)}:] Integer. A two dimensional array of size
  $\mathtt{N\times N_b}$, where $\mathtt{N_b}$ is the number of
  bootstraps, and $\mathtt{N}$ is the range of each bootstrap.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=Bootstrap,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Resampling some data.,
                   label=bootstrap]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 8, Nbt = 5
  Real (kind=DP) :: X(Nmax)
  Integer :: Ib(Nmax, Nbt), I, J


  CALL Random_Number(X)
  ! Generate 5 bootstraps
  CALL BootStrap(Ib)

  ! Write the original sample, and the bootstraps
  Write(*,'(1000ES33.25)')(X(J), J=1, Nmax) 
  Do I = 1, Nbt
     Write(*,'(1000ES33.25)')(X(Ib(J)), J=1, Nmax) 
  End Do

  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{SaveBstrp(Ibt, Filename)}}
\index{SaveBstrp@Subroutine \texttt{SaveBstrp(Ibt, Filename)}}

\subsection{Description}

Saves the bootstrap stored in $\mathtt{Ibt}$ and saves it in the file
\texttt{Filename}.

\subsection{Arguments}

\begin{description}
\item[\texttt{Ibt(:)}:] Integer. A two dimensional array of size
  $\mathtt{N\times N_b}$, where $\mathtt{N_b}$ is the number of
  bootstraps, and $\mathtt{N}$ is the range of each bootstrap.
\item[\texttt{Filename}: ] Character (len=*). A file name to save the
  resampling data.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=SaveBstrp,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Reading the resampling info.,
                   label=SaveBstrp]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 8, Nbt = 5
  Integer :: Ib(Nmax, Nbt)


  ! Generate 5 bootstraps
  CALL BootStrap(Ib)

  ! Save it
  SaveBstrp(Ibt, 'example.bst')

  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{ReadBstrp(Ibt, Filename)}}
\index{ReadBstrp@Subroutine \texttt{ReadBstrp(Ibt, Filename)}}

\subsection{Description}

Reads the bootstrap stored in the file \texttt{Filename}, and returns
it in  $\mathtt{Ibt}$.

\subsection{Arguments}

\begin{description}
\item[\texttt{Ibt(:)}:] Integer. A two dimensional array of size
  $\mathtt{N\times N_b}$, where $\mathtt{N_b}$ is the number of
  bootstraps, and $\mathtt{N}$ is the range of each bootstrap.
\item[\texttt{Filename}: ] Character (len=*). A file name to read the
  resampling data.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=ReadBstrp,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Saving the resampling info.,
                   label=ReadBstrp]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 8, Nbt = 5
  Integer :: Ib(Nmax, Nbt)


  ! Read a saved Bootstrap.
  ReadBstrp(Ibt, 'example.bst')

  Stop
End Program Tests
\end{lstlisting}

\section{Subroutine \texttt{EstBstrp(Data, Ibt, Func, Val, Err[, Rest])}}
\index{EstBstrp@Subroutine \texttt{EstBstrp(Data, Ibt, Func, Val, Err)}}

\subsection{Description}

Estimates using the Bootstrap method the average and error of an
estimator given as a user supplied function.

\subsection{Arguments}

\begin{description}
\item[\texttt{Data(:)}:] Double precision Real. A one dimensional
  array with the original sampling.
\item[\texttt{ibt(:,:): }] Integer two dimensional array. The bootstrap
  that we want to use to make the estimation.
\item[\texttt{Func}: ] A user suplied function that returns the value
  of the estimator. An interface block of the following type should be
  defined. 
\begin{verbatim}
    Interface 
       Function Func(X)
         USE NumTypes
         
         Real (kind=DP), Intent (in) :: X(:)
         Real (kind=DP) :: Func

       End Function Func
    End Interface
\end{verbatim}
\item[\texttt{Val}: ] Double precision real. Output. The value of the 
  estimation of the parameter. 
\item[\texttt{Err}: ] Double precision real. Output. An estimation of
  the error in the estimation of the parameter. 
\item[\texttt{Rest}: ] Double precision real one dimensional array
  (same dimension as the number of bootstraps in
  \texttt{Ibt}). Output. The value of the estimator for each
  resampling. 
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=EstBstrp,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Estimating the average.,
                   label=EstBstrp]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 100, Nbt = 50
  Integer :: Ib(Nmax, Nbt)
  Real (kind=DP) :: Avg, Err, Data(Nmax), Rest(Nbt)

  Interface 
     Function F(X)
       USE NumTypes
         
       Real (kind=DP), Intent (in) :: X(:)
       Real (kind=DP) :: F

     End Function F
  End Interface

  ! Read a saved Bootstrap, and the data from a file
  ReadBstrp(Ibt, 'example.bst')
  Open (Unit=22, File="data.dat")
  Read(22,*)Data
  Close(22)

  ! And estimate the average
  CALL EstBstrp(Data, Ibt, F, Avg, Err, Rest)

  ! Print the Average of each resampling
  Do I = 1, Nbt
    Write(*,*)I, Rest(I)
  End Do

  Stop
End Program Tests

Function F(X)
  USE NumTypes
  USE Statistics
         

  Real (kind=DP), Intent (in) :: X(:)
  Real (kind=DP) :: F

  F = Mean(X)

End Function F

\end{lstlisting}


\section{Subroutine \texttt{BstrpConfInt(Data, Ibt, alpha, Func, dmin, dpls))}}
\index{BstrpConfInt@Subroutine \texttt{BstrpConfInt(Data, Ibt, alpha, Func, dmin, dpls)}}

\subsection{Description}

Gives an Confidence interval for the estimator given as the user
supplied function \texttt{Func}, such that
\begin{displaymath}
  \mathcal P(dmin<Func(Data)<dpls) = 1-2\alpha
\end{displaymath}

\subsection{Arguments}

\begin{description}
\item[\texttt{Data(:)}:] Double precision Real. A one dimensional
  array with the original sampling.
\item[\texttt{ibt(:,:): }] Integer two dimensional array. The bootstrap
  that we want to use to make the estimation.
\item[\texttt{alpha}: ] Double precision real. The level of the
  confidence interval. 
\item[\texttt{Func}: ] A user suplied function that returns the value
  of the estimator. An interface block of the following type should be
  defined. 
\begin{verbatim}
    Interface 
       Function Func(X)
         USE NumTypes
         
         Real (kind=DP), Intent (in) :: X(:)
         Real (kind=DP) :: Func

       End Function Func
    End Interface
\end{verbatim}
\item[\texttt{dmin}: ] Double precision real. Output. The lower limit
  of the confidence interval.
\item[\texttt{dpls}: ] Double precision real. Output. The higher limit
  of the confidence interval.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=BstrpConfInt,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Giving a confidence interval.,
                   label=BstrpConfInt]
Program Tests

  USE NumTypes
  USE Error
  USE Statistics

  Integer, Parameter :: Nmax = 1000, Nb = 10000, Ndiv = 50
  Real (kind=DP) :: Rdata(Nmax), Rmean(Nb), avg, Xmin, Xmax, h, Xac,&
       & err, dmin, dpls
  Integer :: Id(Nmax), Ib(Nb, Nmax), Ntics(Ndiv)

  Interface 
     Function F(X)
       USE NumTypes
         
       Real (kind=DP), Intent (in) :: X(:)
       Real (kind=DP) :: F
       
     End Function F
  End Interface
  
  CALL Normal(Rdata)

  avg = Mean(Rdata)
  ! Now create the resamples
  CALL Bootstrap(Ib)

  CALL EstBstrp(Rdata, Ib, F, avg, Err, Rmean)
  Write(*,*)'#', avg, Err, Mean(Rdata)

  CALL BstrpConfInt(Rdata, Ib, 0.1_DP, F, dmin, dpls)
  Write(*,*)'Intervalo:', dmin, dpls
  Write(*,*)Avg - Dmin, Dpls - Avg
  Write(*,*)(dpls - dmin)/2.0_DP, 1.64485_DP/Sqrt(Real(Nmax,kind=DP))

  Stop

End Program Tests

Function F(X)
  USE NumTypes
  USE Statistics
  
  Real (kind=DP), Intent (in) :: X(:)
  Real (kind=DP) :: F

  F = Mean(X)
  
  Return
End Function F
\end{lstlisting}

\section{Function \texttt{Prop\_Error(X, Dx, Func[, N])}}
\index{Prop@Function \texttt{Prop\_Error(X, Dx, F[, N])}}

\subsection{Description}

Being \texttt{Func} a function of one or several variables, and \texttt{X} a point
known 
with accuracy given by \texttt{Dx}, this routine propagate the errors
in the position of the point, to obtain the accuracy of \texttt{F(X)}.

\subsection{Arguments}

\begin{description}
\item[\texttt{X[(:)]}:] Double (DP) or simple (SP) precision one
  dimensional array or number. The value(s) of the variables.
\item[\texttt{Dx[(:)]}:] Double (DP) or simple (SP) precision one
  dimensional array or number. The value(s) of the error in the
  variables.
\item[\texttt{Func}:] A user supplied function. The function of the
  variables whose errors you want to propagate. An interface block of
  the following type should be defined. 
\begin{verbatim}
    Interface 
       Function Func(X)
         USE NumTypes
         
         Real (kind=DP), Intent (in) :: X(:)
         Real (kind=DP) :: Func

       End Function Func
    End Interface
\end{verbatim}
\item[\texttt{N}:] Integer, Optional. The number of samples used to
  propagate the error. The default value is 1000.
\end{description}

\subsection{Output}

A real double or simple precision (same type as the input). An stimate
of the accuracy of the value \texttt{F(X)}.

\subsection{Examples}

Here we compare our method of propagating errors, with the traditional
propagation of errors formula for the case of the function
$f(x,y)=x^2+y^2$. 

\begin{lstlisting}[emph=Prop_Error,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Propagating errors in a function.,
                   label=properror]
Program Tests

  USE NumTypes
  USE Statistics

  Real (kind=DP) :: Xm(2), Em(2)

  Interface 
     Function F2(X)
       USE NumTypes
       Real (kind=DP), Intent (in) :: X(:)
       Real (kind=DP) :: F2
     End Function F2
  End Interface


  Xm(1) = 1.23_DP
  Xm(2) = 1.62_DP
  Em(1) = 0.02_DP
  Em(2) = 0.03_DP
  Write(*,*)'Value of the function at Xm: ', F2(Xm)
  Write(*,*)'  Error estimate:            ', Prop_Error(Xm, Em, F2)
  Write(*,*)'  Error estimate by hand:    ', &
       & Sqrt((Em(1)/Xm(1))**2 + (Em(2)/Xm(2))**2)*F2(Xm)

  Stop
End Program Tests


Function F2(X)
  USE NumTypes
  
  Real (kind=DP), Intent (in) :: X(:)
  Real (kind=DP) :: F2

  F2 = X(1)**2 + X(2)**2 
       
  Return
End Function F2
\end{lstlisting}

\section{Subroutine \texttt{MCIntegration(X1, X2, Func, Val, Err[, Tol])}}
\index{MCIntegration\texttt{MCIntegration(X1, X2, Func, Val, Err[, Tol])}}

\subsection{Description}

Uses the Monte Carlo method to stimate the value of the integral of
the several variables function \texttt{Func} between the limits
(\texttt{X1(1)}, \texttt{X1(2)}, \dots, \texttt{X1(Ndim)}) and 
(\texttt{X2(1)}, \texttt{X2(2)}, \dots, \texttt{X2(Ndim)}). The value
of the integral is returned in \texttt{Val}, and an estimation of the
error in \texttt{Err}. The desired precision can be introduced via the
optional argument \texttt{Tol}.


\subsection{Arguments}

\begin{description}
\item[\texttt{X1(:)}:] Double or simple precision Real. A one dimensional
  array with the lower limits of the integrals.
\item[\texttt{X2(:)}:] Double or simple precision Real. A one dimensional
  array with the upper limits of the integrals.
\item[\texttt{Func}: ] A user suplied function that returns the value
  of the function to integrate. An interface block of the following
  type should be defined. 
\begin{verbatim}
    Interface 
       Function Func(X)
         USE NumTypes
         
         Real (kind=DP), Intent (in) :: X(:)
         Real (kind=DP) :: Func

       End Function Func
    End Interface
\end{verbatim}
\item[\texttt{Val}: ] Double or simple precision real (same as
  input). Output. An estimation of the value of the integral.
\item[\texttt{Err}: ] Double or simple precision real (same as
  input). Output. An estimation of the error.
\item[\texttt{Tol}: ] Double or simple precision real. Optional. An
  estimation of the desired precision.
\end{description}

\subsection{Examples}

\begin{lstlisting}[emph=MCIntegration,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Integrating a two dimensional function.,
                   label=BstrpConfInt]
MODULE ExtData

  USE NumTypes
  
  Real (kind=DP) :: Sig = 2.345763450_DP
  Real (kind=DP) :: X1 = -0.324563562345643650_DP, X2 = 1.623473_DP, &
       & Y1 = 1.34764574570_DP, Y2 = 5.3476347634_DP


End MODULE ExtData

Program MonteCarlo
  
  USE NumTypes
  USE Error
  USE Statistics

  USE ExtData


  IMPLICIT NONE

  Real (kind=DP) :: Res, Err

  Interface
     Function F(X)
       USE NumTypes
              
       Real (kind=DP), Intent (in) :: X(:)
       Real (kind=DP) :: F
     End Function F
  End Interface

  CALL MCintegration((/X1,Y1/), (/X2,Y2/), F, Res, Err, 1.0E-4_DP)
  Write(*,'(1A,100ES33.25)')'# Value of the integral: ', Res, Err
 

  Stop
End Program Monte


Function F(X)

  USE NumTypes
  USE ExtData

  Real (kind=DP), Intent (in) :: X(:)
  Real (kind=DP) :: F


  F = exp(-1.0_DP/(2.0_DP*Sig**2) * (Sum(X(:)**2)))*Sin(X(1)*X(2))

  Return
End Function F
\end{lstlisting}



% Local Variables: 
% mode: latex
% TeX-master: "lib"
% End: 

